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We consider a cloud of fermionic atoms in an optical lattice described by a Hubbard model with 
an additional linear potential. While homogeneous interacting systems mainly show damped Bloch 
oscillations and heating, a finite cloud behaves differently: It expands symmetrically such that gains 
of potential energy at the top are compensated by losses at the bottom. Interactions stabilize the 
necessary heat currents by inducing gradients of the inverse temperature 1/T, with T < at the 
bottom of the cloud. An analytic solution of hydrodynamic equations shows that the width of the 
cloud increases with t 1 ^ 3 for long times consistent with results from our Boltzmann simulations. 
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Measuring the conductivity is probably the most fun- 
damental experiment when investigating the properties 
of a metal. The influence of constant forces arising from 
electric or gravitational fields on quantum particles in pe- 
riodic potentials is an old and well studied problem [Tj. 
For free particles, Bragg scattering from the periodic po- 
tential induces Bloch oscillations [JJ-[3] . While difficult 
to observe in solids due to disorder and interactions, such 
Bloch oscillations have been measured in semiconductor 
superlattices [J] and for ultracold bosonic atoms in opti- 
cal lattices [2 [3] created by standing waves of lasers, e.g., 
to determine the masses of atoms with high precision [5] . 

Recent experimental developments make it also possi- 
ble to load fermionic ultracold atoms in the lowest band 
of optical lattices. Using an equal population of two hy- 
pcrfine levels and Feshbach resonances (e.g., of 40 K) it 
became possible to realize the Hubbard model (see be- 
low) with tunable interaction strength. While the tem- 
peratures in current experiments are still rather high, 
it was possible to see [HI [7] signatures of the onset of 
a metal-insulator transition induced by strong interac- 
tions. Constant external forces can be realized using ei- 
ther gravitation or accelerated lattices j2 a (after a careful 
elimination of other confining potentials as in Ref. [5]). 

The Hubbard model in an electric or gravitational field 
has attracted previously considerable attention [5lU6|. 
partially motivated by the question how large electric 
fields can lead to a breakdown of the Mott insulating 
state. When discussing the physics of such systems ei- 
ther for weak or strong interactions, it is important to 
take energy conservation into account. While real solids 
are usually probed in contact with some thermal bath, ul- 
tracold atoms provide almost ideal realizations of closed 
quantum systems implying severe restrictions on the dy- 
namics. For a translationally invariant, infinite system 
one expects that even weak interactions lead to a damp- 
ing of the Bloch oscillations (at least in the absence of 
superfluidity and in dimensions larger than 1). During 
this process, potential energy is converted into heat. As 
long as there is no coupling to an external bath which can 




FIG. 1: (color online) Atomic cloud in an optical lattice in 
the presence of gravity. A symmetric expansion of the cloud 
is possible by transporting energy "uphill" as atoms lose po- 
tential energy at the bottom and gain it at the top. This en- 
ergy current is associated with a gradient of 1/T with T < 
(T > 0) at the bottom (top) of the cloud. 

transport the heat away, the system gets hotter and hot- 
ter. It finally reaches a steady state characterized by an 
infinite temperature and vanishing current. For moder- 
ately strong interactions we expect that the same steady 
state is reached, but Bloch oscillations may become over- 
damped. Finally, in the Mott regime, the currents and 
dissipation are exponentially suppressed for small fields 
O [10] , but in the long-time limit heating up to T — oo is 
expected. While we are not aware of a study which care- 
fully tests this plausible picture, it seems to be consistent 
with the available numerical results, e.g. the short-time 
behavior extracted from dynamical mean-field equations 
studied by Eckstein, Oka and Werner [9]. 

Here we argue that the physics of an inhomogeneous 
system, i.e., a finite atomic cloud, is rather different; see 
Fig. [l] For bosons, this question has been studied using 
the Gross-Pitaevskii equation [14] [15] . The motion of the 
system is governed by energy conservation and the fact 
that the kinetic and interaction energies per atom are 
bounded not only from below but also from above (as- 
suming that the optical lattice is sufficiently deep that 
interband transitions can be safely neglected). Energy 
conservation prohibits that the cloud as a whole moves 
up or down over long distances. For noninteracting parti- 
cles this implies that only Bloch oscillations are possible. 
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With interactions a second possibility arises: energy con- 
servation allows that the cloud expands symmetrically, 
losing potential energy at the bottom and gaining it at 
the top of the cloud. To this end, energy has to be trans- 
ported over macroscopic distances upwards. We can ex- 
pect that the interacting system explores this part of the 
phase space: the cloud expands. For a sufficiently large 
cloud and not too strong driving forces, one can expect 
that most of the system is approximately in local equi- 
librium such that a local temperature can be defined. 
The presence of an energy current and the absence of 
a particle current in the center of the cloud implies a 
gradient in temperature or rather in (3 = l/T. Com- 
bining this insight with the expectation that T — > oo in 
the center of the cloud for long times (as in the homoge- 
neous system) , we obtain the situation sketched in Fig. [I] 
The cloud expands, using an energy current driven by 
a gradient of j3 = l/T with T > at the top of the 
cloud and negative absolute temperatures, T < 0, at the 
bottom. The numerical confirmation of this qualitative 
picture by a Boltzmann simulation and the associated 
quantitative theory presented below are the main results 
of this Letter. Note that negative T, i.e. systems where 
high-energy states have a higher occupation, e~ E ' kBT , 
than low-energy states, are well-defined for closed sys- 
tems with an upper bound in energy. In Ref. |17j we 
discuss how T < can be realized and detected. 

To substantiate our picture, we study the dynamics of 
a cloud governed by a two-dimensional (2D) fermionic 
Hubbard model in the presence of a gravitational force, 

H = —J ^2 c i<7 C 3°- + U ^ ^»t n »4- + ^ Vi^itr, (1) 

(ij),<7 i t jCT 

Vt = grf, 

with hopping rate J and interaction U. The number of 
atoms of spin a at site i located at position is given 
by Tii a = c\ a Ci a . We set the lattice constant a = 1 in 
the following. The constant force g is applied in the x 
direction, and this term prohibits equilibrium at finite 
T . As we are mainly interested in the dynamics in this 
direction, we consider for simplicity a translationally in- 
variant system in the y direction. Initially the system is 
in equilibrium at a given T = J (a typical temperature 
for current experiments), confined by an extra harmonic 
trap in the x-direction, which is switched off at t = 0. For 
g = 0, the Hamiltonian has been realized in Ref. [5] 
where the expansion dynamics was studied also theoret- 
ically. Considering a 2D instead of a 3D setup has the 
advantage that in experiments it allows one to change the 
value of g by tilting the 2D lattice. While g sa J can be 
reached easily in experiments jS], we are more interested 
in a regime where g is much weaker. As an alternative 
to tilting, one can also use accelerated lattices 0. 

For not too large values of U and g, the system can be 
described by a semiclassical Boltzmann equation. As for 
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FIG. 2: (color online) Densities n(x, t) (red) and inverse local 
temperatures P(x,t) = 1/T(x,t) (blue dots) from the Boltz- 
mann simulation for U = J and g ~ 0.13 J. The expansion is 
associated with a l/T gradient with T < 0, T — oo, T > at 
the bottom, in the center and at the top of the cloud, respec- 
tively. For long t, f) approaches fi^ 1 ' = — d x n/(nd x V) (green). 
In the tails, local temperatures cease to be a useful concept 
and the noise reflects Bloch oscillations. 



an inhomogeneous system this is a numerically demand- 
ing integro-differcntial equation, we employ a version of 
the relaxation time approximation [8 where both particle 
number and energy is conserved: 

c\/ + v k .V r / + F.V k / = --^-(.f-/° e ), (2) 

r(n,e) 

Here /(r, k, t) is the occupation probability in phase 
space and the force term F = — V r (g x) — lTV r n(r, t) con- 
tains the external potential and interaction corrections on 
the Hartree level. The velocity is defined as v k = Vk£ k 
with the dispersion relation e k = — 2 J[cos(k x ) +cos(k y )]. 
The parameters T(r, t) and /i(r, t) of the Fermi function 
/° e are determined from n = J d 2 k /°(e k — /j,, T) and 
e = 4^2 / d 2 ke k f°(e k - (i,T), where n = n t = n i = 
n(r,t) and e = e(r, t) are the local particle and kinetic 
energy densities per spin, respectively. We use a scatter- 
ing rate l/r(n, e) ~ U 2 which arises from interparticle 
scattering with momentum transfer to the lattice (umk- 
lapp) and which was determined in Ref. [5j to reproduce 
the e- and n-dependent diffusion constant of the Hub- 
bard model to order U 2 (obtained from an independent 
calculation). In the regime of high T, relevant for our 
study, l/r(n, e) w n(l — ti)/tq is approximately indepen- 
dent of temperature with 1/tq « 0.609 U 2 / J in D = 2 
(see Ref. [5] for more details). 

A typical result of the Boltzmann equation is shown 
in Fig. [2] In the gravitational potential, the cloud falls 
initially until energy- and temperature profiles become 
antisymmetric relative to the center of the cloud such 
that T = oo in the center. For the parameters of Fig. [2] 
Bloch oscillations of the center xo(t) ~ Jdx xn are over- 
damped but become visible for weaker interactions; see 
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inset of Fig. [3] Note that there are always Bloch oscilla- 
tions in the tails of the cloud where low densities prohibit 
scattering. The most prominent feature is, however, the 
symmetric expansion of the cloud in the long-time limit. 
The width of the cloud increases slowly; see Fig. [3j For 
the local temperatures (defined above) one observes the 
formation of a gradient in 1/T with negative (positive) 
absolute temperatures at the bottom (top) of the cloud. 
As discussed above, this gradient is associated to the heat 
currents needed for the symmetric expansion of the cloud. 

To obtain an analytic theory of the effects described 
above, we derive the hydrodynamic equations for n(r,t), 
and e(r, t) in the regime where the scattering time is 
smaller than the Bloch oscillation period, gr /2ir <C 1. 
The opposite limit of weakly damped Bloch oscillations is 
briefly discussed in the conclusions. We use the standard 
approach to expand the Boltzmann Eq. ([2| to linear order 
around the local equilibrium distribution f~ n (ek). As we 
are interested in the dynamics for long times, we use that 
T will become very high and n small with 



3000 



T = ~AJ 2 n/e 

to leading order. Using also that 1/t ; 
we obtain (see supplement |18j ) 
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where j„ and j e are particle and energy currents. Above 
we have included all subleading terms in the high-T 
expansion, which contribute in the t — > oo limit ac- 
cording to the analysis below. Kinetic energy is cre- 
ated by Joule heating with the rate — jnVU(x) where 
V(x) = gx + Un(x) includes Hartree corrections. 

Note that the hydrodynamic equations can be also de- 
rived directly from the Hubbard model. Up to changes 
in numerical prefactors, Eqs. Q and ([5| therefore should 
be exact for the large T limit of the Hubbard model in 
dimensions D > 1 for large clouds and weak potentials 
V (the breakdown of the hydrodynamic equations found 
for g — in [5] is not important here). In that sense, 
our analysis and results below do not rely on our specific 
Boltzmann equation. Also note that the ID case is more 
subtle due to the integrability of the Hubbard model. 

Surprisingly, it is possible to analyze the long-time 
limit of the nonlinear, coupled partial differential Eqs. Q 
and ([5| analytically. Measuring distance x relative to the 
center off mass, we start from the scaling ansatz 



n(x,t)=N —F[x/R(t)} 



(6) 



where F is a dimensionless scaling function of unit width 
and Jdz F[z] = 1 such that Nq — Jdx n, R 2 (t) — 
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FIG. 3: (color online) To a good approximation, the cloud 
radius cubed grows linearly in time, R 3 (t) ~ (t — to), for a 
wide range of U and g. Inset: Center of mass of the cloud 
Xo(t) showing damped Bloch oscillations for weak U. 



Jdx x 2 n/No. Our goal is to calculate both R(t) and F[z] 
in the long-time limit where R(t) is large. On the top 
panel of Fig. [3] we show the scaling function F obtained 
by rescaling the density profiles from the Boltzmann sim- 
ulation according to Eq. jfjj. To obtain a scaling ansatz 
for e, the most important step is to realize that to lead- 
ing order, energy conservation prohibits expansion. We 
therefore set e = eo + Ae, and determine eo by setting 
j% — in Eq. ^ for e = eo (neglecting subleading terms 
cx e 2 ). This leads to the ansatz 



eo 



U 2 



Or 



eo + ^i+^ G l x / R (^ ( 7 ) 



with an unknown exponent 7 > 1 (such that Ae <C eo) 
and a new scaling function G. To check whether e eo to 
leading order approximation, we compare 1/T from our 
numerical data to Eq. ([3]) with e ~ e and find excellent 
agreement (see Fig. [2|. While this provides already a 
precise, quantitative theory for the temperature gradient, 
only corrections to eo determine how the cloud expands. 

From the particle number continuity equation in 
Eq. (jij) we can calculate j% = — J x d t n(x' ,t)dx' = 
N xFR~ 2 R using Eq. (jH]), which has to be equivalent 
to 2%, m Eq. ([5| . We can match the scaling dimension of 
this term to the subleading terms proportional to Ae if 
R~ f ~ 1 R = const., implying that 



R(t) 



a 



\N g 2 



1/7 



(t-t ) lh 



(8) 



where ao and to are integration constants and the other 
prefactors of R(t) have been chosen for later convenience. 

We can now substitute these results into the energy 
continuity equation, Eq. Q, from which we obtain terms 
proportional to t' 1 ,t~ 3 ^ ,t~ 1 - 2 ^ and t-' 2 - 1 ^. From 
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FIG. 4: (color online) Top: Normalized density as a function 
of z = x/R(t) for tJ = 500, 1500, 2500, 3500. As predicted by 
Eq. (10 1 the shape is approximately Gaussian. Middle: En- 
ergy density e compared to eo given by Eq. Q. Bottom: 
Comparison of the subleading correction, Ae, to Eq. ( |12[ ) 
where j" and n are taken from our simulation. 



period, r <C is = 27r/g, we derived analytically a precise 
quantitative theory for the expansion of the cloud based 
on hydrodynamics. It is an interesting question for fur- 
ther studies to investigate analytically also the opposite 
limit, r 3> ts, where diffusion constants are proportional 
to 1/r instead of r as only scattering allows the atoms 
to move forward. Nevertheless, energy conservation will 
again lead to a symmetric expansion, however, with re- 
duced speed. For t — > oo one reaches this regime for 
arbitrary U as for t ^> Nq /(a^gJ A TQ ) one always obtains 
t 3> ts- For moderately strong interactions, however, 
we expect that atomic losses will make it very difficult to 
track experimentally the system for such long times. 

We acknowledge financial support by SFB TR 12 and 
SFB 608 of the DFG and the German National Academic 
Foundation and useful discussions with I. Bloch, E. Dern- 
ier, V. Gurarie, M. Kunze, S. Manmana, A. M. Rey, U. 
Schneider, and R. Schutzhold. 



the condition that the two leading terms have to cancel 
to fulfill energy conservation, we obtain 

3, R~t 1/3 
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Therefore, we conclude that the width of the cloud in- 
creases with i 1 / 3 for long times in the diffusive regime. 
In comparison, subdiffusive exponents 7 -1 < 1/4 |15j 
and 7 _1 = 0.19 [14] have also been obtained for bosons 
within time-dependent mean-field theory. The parame- 
ter «o and therefore the shape of F depend on the initial 
conditions but for our simulations we find that «o/18 is 
small and therefore the rescaled density profile F is close 
to a Gaussian according to Eq. (10 1. 



Figure [4] shows that the analytic approach describes 
the results of the Boltzmann simulations with high pre- 
cision. We combine Eqs. ^ and (11) to predict 
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g3 n 2 



g3 n 2 



(12) 



Testing this relation numerically, we find good agreement 
even for this subleading term; see Fig. [4] As expected, 
our hydrodynamic equations are, however, not valid in 
the tails of the cloud where 1/r is small. 

In conclusion, we have shown that interactions lead to 
a symmetric expansion of fermionic clouds in optical lat- 
tices subject to a constant force. In the limit where the 
scattering time is short compared to a Bloch oscillation 
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THE MOTION OF THE CENTER OF MASS OF 
THE CLOUD 

The center of mass of the cloud is defined by 

x (t) = —r [dx xn(x,t) , (13) 
J 

where N = Jdx n(x,t). Using energy conservation, 

dx ^e(x, t) + gxn(x, t) + -Un 2 (x, t)j — const. (14) 

we can estimate Xo(t) for very short and very long times. 
With xo(t = 0) = 0, the center of mass of the cloud at 
time t is simply given by 



x o(t) = ^t— dx(e(x,0) - e(x,t)) 



U 



2N Q g 



dx(n 2 (x,0) - n 2 (x,t)) . (15) 



As is shown in Fig. [5] for short times, the total Hartree 
energy decays only very slowly (see below) . For short (or 
moderately long) times, we can therefore approximate 
it as a constant. For the case of weakly damped Bloch 
oscillations it is useful to examine Eq. |l5]) for t = n/g, 
i.e., after half of a Bloch period. In the non-interacting 
case, for all particles k x (t = n/g) = k x (t = 0) + n while 
k y (t = n/g) = k y (t = 0). In d = 2 this implies e{x,t — 
n/g) — as the two terms in the kinetic energy cancel 
each other. As a consequence, 



Xo{t = n/g) « x 



1 

Nog 



dx e(x, 0) 



(16) 



also for weak interactions. For high initial temperatures 
T(t = 0) > J we can approximate e « —AJ 2 n/T and 
obtain x w -4J 2 /(Tg). For T(t = 0) = J used in our 
simulations, one gets Xq ~ ~ 1-59 J/g, which describes the 
numerical results to good approximation, see the inset in 
Fig. 3. of the main text. 

The same arguments can be applied to t n — (2n + 
l)n/g, n £ Z for weakly damped Bloch oscillations. This 
explains the form of Xo(t) shown in the inset of Fig. 3 
(main text) for small U and/or large g, see also Fig. [5] 
While the amplitude of the Bloch oscillations is damped, 
the local minima of Xo(t), are at an approximately con- 
stant value, Xo(t n ) sa x . Therefore, on the time scales 



where on the one hand the Bloch oscillations have died off 
and on the other hand the change of the Hartree energy 
can still be neglected, we obtain 

x (t) rj x . (17) 

Formally, the Hartree energy vanishes for t — > oo and 
therefore 



Xo(t — > oo) = Xq + 



U 



2N g 



dxn 2 (x,0). (18) 



This asymptotic value is, however, reached only very 
slowly as Jdxn 2 cx 1/R(t) oc 1/i 1 / 3 in the regime where 
our hydrodynamic approach is valid. 

For the validity of the scaling ansatz used in the main 
text, where an antisymmetric form of e(x, t) = — e(— x, t) 
was used, it is important to check that the total kinetic 
energy vanishes much faster than the Hartree energy 
which is indeed the case, see Fig. [5] 
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FIG. 5: (color online) Total kinetic energy, potential energy, 
and Hartree energy (per particle) for different times as a func- 
tion of t for U/J = 1, g/J = 0.325 and g/J = 1.3. 



DERIVATION OF HYDRODYNAMIC 
EQUATIONS FROM THE BOLTZMANN 
EQUATION 

In this section we shall derive hydrodynamic equations 
from the Boltzmann equation 

a t / + v k -V r /-(V r ^)-V k / = -r- 1 (n,e)(/-/ ) (19) 



() 



for weak potentials V and large scattering rates r _1 for 
high temperatures. The local densities and currents per 
spin component are defined by 

n = i / d2k f > e = 4^J d2k £k/ ' (20) 

Jn = Je = i/ rf2k V k e k /.(21) 

We start from the particle number and kinetic energy 
continuity equations, 



n = -V r j„ , 

e - -V r j e + (-V r F)-j„ 



(22) 



Note that although one should also consider the term 
(— r)dt/ [ra(r, f),e(r, t), ek] in Eq. (23), it does not con- 
tribute to the currents due to the odd power of v k in the 



corresponding momentum integral in Eq. (21 1 



For high temperatures f3 J <C 1, the local kinetic energy 
density e is small compared to the bandwidth. Thus we 
expand the distribution function in j3 J, 



1 



c 



(-1+0C 2 . 2 

i + C (i + C) 2 ^ ' 2(1 + C) 3 kP 



(26) 



Our goal is to express the currents j„ and j e as functions 
of n and e in lowest orders in e. 

In the diffusive limit, where rg -C 1, we can solve the 
Boltzmann equation iteratively to get 

/ « f° + (-r)(v k • V r /° + (-V r V) ■ V k /°) , (23) 
where 



1 



C{r,t)e^ r ^ + 1 



(24) 



with the fugacity £ and inverse temperature (3 chosen 
such that 



1 

47r2 



Jd 2 k /°, e=~ Jd 2 k e k /°. (25) 



which we substitute in the definitions of n and e. We can 
easily perform the momentum integrals, and find 



i + C 
-A 2 (d) 



A 2 (d) 



(-1 + 0C .2 
2(1 + C) 3 ' 



(27) 



(i + O 



where A„(<2) = Jd 2 k e™, e.g., A 2 (cZ = 2) = 4J 2 and 
Ai(d = 2) = 36 J 4 . We can solve these equations for £ 
and /3 as a function of rt up to order e 3 . Substituting 
these results back to f° gives 



ee k 

A a (d) 



e 2 (l - 2n)(A 2 (d) - e k ) e 3 (l - 6(1 - n)n)e k (e k " A 4 (d)/A 2 (d)) 



2A|(d)(l -n)n 



+ 



6A 2 (d) 3 (l-n) 2 n 2 



(28) 



We use this expression in Eq. (23) to get / = f[n, e]. Then from Eq. (21 ) we obtain the currents by keeping terms up 
toO(e 2 ): 



-t C V r n + r 



C 



(V r 10 e 



A 2 (d) 

A 2 {d)C - B 2 (d) (1 - 2n+ 2n 2 ) 



2^1 (d) 

B 2 (d) 



A 2 (d) 



V r e + t 



(1 — n) 2 n 2 
B 2 (d) l-2n 
(l-n)n 



(V r n)e 2 



(VrV)£ 



_A 2 (d)C-B 2 (d) (l-2n) 



(1 — n)n 



eV r e 



(29) 
(30) 



Here we used that ^ /d 2 k «k u k = %C; C = 2J 2 a 2 and 
4^2 JcPk = SijB 2 {d);B 2 (d = 2) = 6J 4 a 2 . Note 

that the first term in Eq. p3| does not contribute to the 
currents since f° corresponds to an equilibrium distribu- 



tion function. Thus the currents are always proportional 
to t w To/n. 

For the hydrodynamic equations in the main text, we 
only keep terms which contribute for t oo. 



